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Abstract 

We have implemented a parallel multigrid solver, to solve the initial data problem for 3 + 

1 General Relativity. This involves solution of elliptic equations derived from the Hamiltonian 

(•~^ ■ and the momentum constraints. We use the conformal transverse-traceless method of York and 

^S] ! collaborators [l|, y, y, |j, |5| which consists of a conformal decomposition with a scalar (f) that 

Q^' adjusts the metric, and a vector potential w'' that adjusts the longitudinal components of the 

<^' 

extrinsic curvature. The constraint equations are then solved for these quantities (j), w'' such 

in ; 

^^ . that the complete solution fully satisfies the constraints. We apply this technique to compare 



o 

u 



X 



with theoretical expectations for the spin-orientation- and separation-dependence in the case of 



> 

f— ^ ' spinning interacting (but not orbiting) black holes. We write out a formula for the effect of the 



Tjlj- . spin-spin interaction which includes a result of Wald |y| as well as additional effect due to the 

o. 

\^ , rotation of the mass quadrupole moment of a spinning black hole. A subset of these spin-spin 



effects are confirmed via our numerical calculations, however due to computer time limitations the 
full parameter space has not yet been surveyed and confirmed. In particular, at the relatively small 
separations (d < 18m) we are able to consider, we are unable to confirm the expected asymptotic 
fall-off of d~^ for these effects. 
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1. INTRODUCTION 

Simulation of binary black hole mergers will play an important part in the prediction, 
detection, and the analysis of signals in gravitational wave detectors. In the usual approach 
to computing the merger of black holes (generically called the 3 + 1 method), one has first to 
initiate the simulation by producing consistent data. Four of the components of the Einstein 
equation do not contain time derivatives of the spatial metric, nor of the momentum of the 
3-metric. These components. Goo = and Gjo = 0, are thus called constraint equations, 
and they must be satisfied in any specification of initial data. (We are interested in black 
hole interactions, which are vacuum, i.e. matter-free, so the right side of the Einstein 
equation is zero: G,,i, = 0.) As we recall in the next section, a conformal decomposition 
fl fl fl fl B B B an„w.l solution of tKe.e components to be put in t.e f„™ of a 
set of four coupled elliptic equations. These elliptic equations are the subject of our work. 
We solve them via a multigrid method which applies concepts from J9| to this problem. 
We demonstrate the accuracy of our data by considering features discussed analytically by 
Wald p]. Wald described the spin-spin effect on the binding energy of two black holes in 
an analytic perturbation scheme, where one hole is much more massive than the other. Our 
computational technology is well suited to simulating these effects for equal mass black holes, 
and we demonstrate agreement in some aspects of the computational spin-spin interactions 
with the analytic estimate, for separations that are not small. 

2. 3 + 1 FORMULATION OF EINSTEIN EQUATIONS 

We take a Cauchy formulation (3+1) of the ADM type, after Arnowitt, Deser, and 
Misner |lO|. In such a method the 3-metric gij and its momentum Kij are specified at one 
initial time on a spacelike hypersurface, and evolved into the future. The ADM metric is 

ds^ = -{a^- (3i(3') dt^ + 2(3i dt dx' + gij dx' dx^ (1) 



where a is the lapse function and /3* is the shift 3-vector; these gauge functions encode the 
coordinatization.f 

The Einstein field equations contain both hyperbolic evolution equations and elliptic 
constraint equations. The constraint equations for vacuum in the ADM decomposition are: 



H=h,R-K,,K'^^K''\ 



0, 



(2) 



H' = Vj {K'^ - g'^K) 



0. 



(3) 



Eq. (J2|) is known as the Hamiltonian constraint; Eq. is the momentum constraint (three 
components). Here R is the 3-d Ricci scalar constructed from the 3- metric, and Vj is the 
torsion- free 3-d covariant derivative compatible with Qij. Initial data must satisfy these 
constraint equations; one may not freely specify all components oi Qij and Kij. 
One of the evolution equations from the Einstein system is 



gij = -2aKij + V j(3i + \/i/3j, 
and this will prove useful in our data setting procedure below. 



(4) 
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le initial value problem have been addressed in the past by several groups 
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It is the case that until 



recently, most data have been constructed assuming that the initial 3-space is conformally 
flat. The method most commonly used is the approach of Bowen and York [5|, which chooses 
maximal spatial hypersurfaces (for which the quantity K = K'^a = 0), as well as taking the 
spatial 3-metric to be conformally flat. 

The chief advantage of the maximal spatial hypersurface approach is numerical simplicity, 
as this choice decouples the Hamiltonian constraint from the momentum constraint equa- 
tions. Besides, for /C = 0, if the conformal background is flat Euclidean 3-space, then there 
are known Kij that analytically solve the momentum constraint [5|. The constraints then 



[I] Latin indices run 1, 2, 3 and are lowered and raised by gij and its 3-d inverse g^^ . The time derivative will 
be denoted by an overdot (') 



reduce to one elliptic equation for the conformal factor 0. Very recently substantial success 



2J|. However, 



has been achieved evolving Bowen-York data using "puncture" methods j2c 

we generally use an alternative choice of background 3-nietric, which is based on a metric 

constructed from single black hole Kerr Schild data 22] ; multiple black holes are constructed 



by a superposition in the conformal background. It has been shown that this process, while 



not exact for multiple black hole data, does contain much of t 



for a single black hole, even a spinning or boosted black hole 25| 



le physics. It clearly is exact 



3.1. Kerr Schild Black Holes 



The Kerr-Schild 2j| form of a black hole solution describes the spacetime of a single 
black hole with mass, m, and specific angular momentum, a = j/m, in a coordinate system 
that is well behaved at the black hole horizon: 



ds^ = r]^^ dx^ dx'' + 2H{xy^l^ dx" dx", 



(5) 



where 77^,^ is the metric of fiat space, if is a scalar function of x'^, and /^ is an (ingoing) null 
vector, null with respect to both the fiat metric and the full metric. 



v^'lf.lu = g^'^i^lu = 0. 



(6) 



Comparing the Kerr-Schild metric with the ADM decomposition Eq. (0), we find that 
the t = constant 3-space metric is: Qij = 6ij + 2Hlilj. Further, by comparison to the ADM 
form, we have 



A = 2i//o/j 



and 



Explicit 



a 



y/l + 2HJi 



(7) 



(8) 



brms of H{x^) and la{x'^) for Kerr black holes are given in a number of references. 
See J2^, 25[, 35[- Many details of the algebraic manipulation of the Kerr-Schild form are 



found in reference 

The extrinsic curvature can be computed from Eq.(@]): 



K, 



1 

2^ 



[Vj(3i + V,/?j - Qij] 



(9) 



Each term on the right hand side of this equation is known analytically; in particular, for a 
black hole at rest, gij = 0. 

3.2. Boosted Kerr-Schild black holes 

The Kerr-Schild metric is form-invariant under a boost, making it an ideal metric to 
describe moving black holes. A constant Lorentz transformation (the boost velocity, v, is 
specified with respect to the background Minkowski spacetime) A"^ leaves the 4-metric in 
Kerr-Schild form, with H and l^ transformed in the usual manner: 



X 
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A^x", (10) 

H'ix"') = H ((A'l)"^ x'^) , (11) 

I'six"^) = A^s I, {{A-'r, x'^) . (12) 

Note that /[, is no longer unity. As the initial solution is stationary, the only time dependence 
comes in the motion of the center, and the full metric is stationary with a Killing vector 
reflecting the boost velocity. The boosted Kerr-Schild data exactly represent a spinning 
and/or moving single black hole. 

3.3. Background data for multiple black holes 

The structure of the Kerr-Schild metric suggests a natural extension to generate the 
background data for multiple black hole spacetimes. We first choose mass and angular 
momentum parameters for each hole, and compute the respective H and /" in the appropriate 
rest frame. These quantities are then boosted in the desired direction and offset to the chosen 
position in the computational frame. The computational grid is the center of momentum 
frame for the two holes, making the velocity of the second hole a function of the two masses 
and the velocity of the first hole. We compute the individual metrics and extrinsic curvatures 
in the coordinate system of the computational domain: 

Adij = Vij + 2 aH aU aIj, (13) 

AKr = ^ Ag""' (V, A(3^ + V, aPj - Ag^,) ■ (14) 



The pre-index A labels the black holes. Background data for A^ holes are then constructed 
in superposition: 

N 

9ij = Vij + X] ^ aHaIi aIj, (15) 

A 

N 

K = J2 ^^^'' (16) 

A 

^ / 1 \ 

A tilde ( ~ ) indicates a background field tensor. Notice that we do not use the attenuation 
functions introduced by Bonning et al.J2y]. 

To give the reader a feel for how closely the Kerr-Schild superposition data resemble a true 
binary black hole spacetime, in Figure^we provide a graph comparing the superposed Kerr- 
Schild background data with the subsequent solutions of the constraint equations (described 
below) . 

4. GENERATING THE PHYSICAL SPACETIME 

We will consider in this paper physical applications which use superposed Kerr-Schild 
backgrounds. When multiple black holes are present, the background superposed Kerr- 
Schild data described in the previous section are not solutions of the constraints, Eqs. 1^- 
(j21). Hence they do not constitute a physically consistent data set. A physical spacetime 
can be constructed by modifying the background fields with new functions such that the 
constraints are satisfied. We adopt the conformal transverse-traceless method of York and 
collaborators [l| which consists of a conformal decomposition and a vector potential that 
adjusts the longitudinal components of the extrinsic curvature. The constraint equations 
are then solved for these new quantities such that the complete solution fully satisfies the 
constraints. We do not consider tiK = 0, nor conformally fiat, solutions. 

The physical metric, Qij, and the trace- free part of the extrinsic curvature, Aij, are related 
to the background fields through a conformal factor 

9ij = (f)'^9ij, (18) 

A'^ = 0-iO(i^J + (z^)'^)^ (19) 
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FIG. 1: Comparison of background Kerr-Schild superposition data (dashed lines) with the final 
output of our elliptic constraint equation solver (solid lines) . We see that the background is quite 
close to the physical solution. These particular data were generated for two holes located at 
X = ibSm, with spins ai = 02 = 0.5, with the spin of the x = —5m hole tipped by rotation about 
the X axis by 9i = Tvr/S, and an excision radius of 0.75m. 

where is the conformal factor, and (Iw) will be used to cancel any possible longitudinal 
contribution to the superposed background extrinsic curvature, w"^ is a vector potential, and 



,~ x*J 






(20) 



The trace K is taken to be a given function 



K = K. 
7 



(21) 



Writing the Hamiltonian and momentum constraint equations in terms of the quantities in 
Eqs. ffTH )) - (PT|) . we obtain four coupled eUiptic equations for the fields and w* [l|: 

V20 = (l/8)(i?0+^A^V- 

r'iA'^ + ilwy^)iA, + (Iw),,)), (22) 

VjilwY^ = '^g'^4>^VjK-VjA'^. (23) 

Our outer boundary condition for 0, namely 

9p(p(0-l))|p^oo = O. (24) 

enforces — * 1 at oo, but does not specify the size (or sign) of the p~^ term in 0. (Here 
p"^ = x^ + y"^ + z^ .) We also take as boundary conditions for the vector w"": 

dpifrw'n,) = 0, (25) 

d,{p^w\S,,-n,n,)) =0, (26) 



where rii is the outward pointing unit spatial normal. Condition (jzlj) is a Robin condition 
commonly used for computational conformal factor determination. Conditions ()25|) and ()26|) 
were derived by Bonning et al.(2fil|. 

5. NUMERICAL METHODS 

We first discuss the computational code and tests, and some code limitations. 

The constraint equations Eq. (j22j) . Eq. (j23j) are solved with a multigrid solver J9|. The 
present code is essentially the same as that described in J9|, except that it has been extended 
to the full set of constraint equations, non-flat backgrounds, and features parallel processing. 
The multigrid scheme is essentially a clever means of eliminating successive wavelength- 
components of the error via the use of relaxation at multiple spatial scales. It makes use 
of some sort of local averaging procedure (e.g. Gauss-Seidel relaxation). Such relaxation 
is extremely effective at eliminating short-wavelength components of the error, or in other 
words, at "smoothing" the error (i.e., the residual, see below). However, relaxation fails 
to operate efficiently on long-wavelength components of the error (components that involve 
discretization points more than a few away from the point at which the solution is sought). 

8 



Multigrid addresses the solution repeatedly on grids of different discretization, achieving the 
same efficiency at smoothing every scale. 

Because the implementation in the method is also in [9|, we do not repeat it here. 

6. MULTIGRID WITH EXCISED REGIONS 

In our formulation, the black holes are represented by excised regions. Because we work 
in Cartesian coordinates, and because we want completely general implementation, we do 
not typically expect that the excision will be defined by overlapping points on the various 
grids of different resolution. 

Our definition of the excision region is that on each grid, the inner boundary consists 
of points that lie just inside, i.e. up to one grid point inside, the analytic location of the 
inner boundary, as shown in Figure El While there are exceptional configurations such as 
cubic excision defined so that the excision boundary lies on points of the coarsest grid, this 
definition means that generically the size of the excision is larger on the finer grids. 

This definition of the inner boundary affects the way in which data are restricted from fine 
grids to coarse grids. Away from the inner boundary, weighted restriction is performed, as 
shown in left pane of Figure El However, if any of the points used in the weighted average lie 
on an inner boundary, then these points are not used and instead a simple "copy" operation 
is performed as shown in the right pane of Figure El The inner boundary points themselves 
may need to be filled in on coarse grids (since on the fine grid they may be excised), and 
to do this we apply a weighted "inward extrapolation" using a parabolic fit to surrounding 
fine grid points. 

This scheme has been implemented in a parallel computing environment, using MPI 
to communicate between processors. Each processor handles a part (a patch) of the total 
domain. The patch is also logically surrounded by "ghost zones". Because we deal with 
a finite-difference representation of derivatives, the communication between processors re- 
quires the filling of these "ghost zones" on the borders of the patches, using values computed 
on other processors, so that derivatives can be accurately computed near the boundaries of 
the patches. This has implications for the way that smoothing is handled in our simulation. 

On a single processor Gauss-Seidel smoothing proceeds across the grid, and the updates at 
any particular point involve some surrounding points that have been updated and some that 
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FIG. 2: Example of how the inner boundary is defined, showing points on a coarse grid and a fine 
grid. Inner boundary points are those points which are immediately interior to a circle of radius rex- 
The large filled circles show normal interior grid points (i.e., non-excised, non-boundary points) 
on the coarse grid, and the large open circles show boundary points on the coarse grid. The small 
filled and open circles show fine grid interior points and boundary points, respectively. The small 
dots show excised points on the fine and/or coarse grids, as appropriate. (Only one quadrant of a 
full domain is shown in this picture for purposes of clarity.) 
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FIG. 3: 1-D schematic of scheme for restriction scheme near inner boundary. The circles on the 
rightmost X's indicate that this is where the Dirichlet conditions are applied, i.e., there are data 
on these points. One could use these points in weighted restriction even in the case shown in the 
right panel. However, we choose never to use these boundary data in weighted restriction, and 
instead do a simple "copy" operation. The boundary points themselves are updated using either 
a direct copy, or for coarse grid points over excised fine grid points via an average over parabolic 
fits of fine-grid data in all available directions. 



10 



have not been. If the same scheme has been implemented on two processors (say sphtting the 
X— axis), the buffer region of one patch will already have been updated when the smoother 
of the other patch begins to use the equivalent points. The order and direction of the filling 
of the ghost zones can lead to inconsistent behavior {i.e., the result will be different from the 
single processor result). One solution to this is to insert "wait" commands into the parallel 
code, so that processors wait to carry out the process in the correct order. This has the effect 
of slowing the execution, and loses the advantage of parallel processing. A better approach 
is to use something like red-black Gauss-Seidel. (In 2-d the red-black pattern is like that on 
a checkerboard.) If the differential operator involves only diagonal second derivatives (no 
mixed partials) then each point is updated using only points of the opposite color. Then 
all the reds can be updated before any of the blacks, and vice versa. This ameliorates the 
ghost zone synchronization problem; the ghost zones can be maintained in the correct state 
for every step. In this case parallelization works as anticipated. 

If the background is taken as flat space, then these conditions apply. But we work with 
Kerr-Schild forms of the metric which guarantee that there will be mixed partial derivatives 
in the operators, and the parallel synchronization problem reappears. 

Our solution is to introduce what we call rainbow smoothing, in which we make a total 
of eight passes (like the two passes in red-black smoothing) over the grid, where each pass 
has a stride width of two over each of the three dimensions of the grid. 

7. VERIFICATION OF CONSTRAINT SOLUTION 

To verify the solution of the discrete equations, we have examined the code's convergence 
in some detail. We use a set of completely independent "residual evaluators" for the full 
Einstein system (here applied only to the initial data), originally constructed by Anderson 
5ll | . These evaluate the Einstein tensor, working just from the metric produced by the com- 
putational solution, to return fourth order accurate results. They are completely different 
from the way the equations are expressed in the constraint solver code. 

Figure m shows such a plot of convergence for the Hamiltonian constraint in an equal mass 
binary black hole spacetime. The holes are located at ±5m on the x-axis, where m is the mass 
of one hole. The elliptic equations were then solved on grids of sizes 385^, 449^, and 513^, 
giving finest-grid resolutions of approximately m/12.8, m/15, and m/17. We use a five-level 
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FIG. 4: Convergence of the Hamiltonian constraint along the positive x and y axes. We show 
values of the constraint obtained and three different finest-grid resolutions, 385^, 449'^ and 513^, 
scaled by appropriate ratios of the mesh spacings consistent with second-order convergence. We 
see that there is good convergence everywhere except near the outer boundaries. Because of this 
loss of convergence near the outer boundaries, we evaluate the ADM mass over the surface a cube 
with half-width 12M. (In the left pane, the vertical scale has been exaggerated in order to zoom in 
on the "body" of the domain, and cuts out the peaks immediately adjacent to the excision regions. 

multigrid hierarchy. Figure |3] demonstrates almost perfect second order convergence, except 
near the outer boundary, where the convergence is apparently first order. The second order 
convergence shows that we have achieved a correct finite difference solution to the initial 
data problem. 

8. COMPUTATIONAL LIMITATIONS ON GRID COARSENESS IN EXCISED 
BLACK HOLE SPACETIMES 

In the examples given here, where we work on a fixed spatial domain (say ±15m), the 
finest resolutions of 385^, 449^, 513^, used in the convergence test correspond to coarsest-grid 
resolutions of 25^, 29^, 33^, respectively. For the 33^ grid, the ±15m domain is discretized 
at about Im resolution. 

The problem of required resolution for black hole simulations has been discussed in this 
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context at least since the early Grand Challenge 28J efforts. Present computational resources 
allow much larger grid size than in the Grand Challenge epoch, so the conffict appears at 
higher resolutions and larger physical domains than previously, and we can do substantial 
physics with the present configuration. Our approach will be to introduce a multiresolution 
scheme to maintain required resolution near the central "action", and allow coarsening 
further away. To accomplish this we are investigating a mesh-refined multigrid, similar to 



that described by Brown and Lowe|29j. However, for the present work, we simply use very 
high resolution, the highest that we can presently achieve on the computers available to us, 
namely 513^ points using 32 processors. 

9. SPIN-SPIN EFFECTS IN BLACK HOLE INTERACTION 

Wald p] directly computes the force for stationary sources with arbitrarily oriented spins. 
He considered a small black hole as a perturbation in the field of a large hole. The result 
found for the spin-spin contribution to the binding energy is 

Here, S, S' are the spin vectors of the sources and fi is the unit vector connecting the 
two sources, and d is any reasonable measure of separation that approaches the Euclidean 
distance x(l + 0((i~^)) at large d (such as the distance measured in the flat background 

using a definition of intrinsic mass that differs 



used in the initial data setting). Dain 

from ours (see below), finds binding energy which agrees with Wald's (Eq. ()27|)) at 0{d~ 
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9.1. Computational Spin-Spin Effects in Black Hole Binding Energy 

In order to investigate a computational implementation validating Eq. (j?7|) . we begin 
with a standard definition of the binding energy for black hole interactions. 

The total gravitational energy in a binary system can be computed from the initial data 
using the ADM mass Madm, which is evaluated by a distant surface integral (see Eq. ()29|1 
below), and gives the Newtonian gravitational mass as measured "at infinity" . For a measure 
of each hole's intrinsic mass, we use the horizon mass Mah defined by (Eq. ()30|)) below. Thus 
the binding energy, Eb, is defined as 

13 



Eb = Madm - Mah - M^H- (28) 

The ADM mass is evaluated in an asymptotically flat region surrounding the system of 
interest, and in Cartesian coordinates is given by 

1 / fdgji dgjj 



The apparent horizon is the only structure available to measure the intrinsic mass of a 
black hole. I Complicating this issue is the intrinsic spin of the black hole; the relation is 
between horizon area and irreducible mass: 

Ah = lQiim?i„ = Sirm (m + \/ [m^ — a^) ) . (30) 

As Eq. ()30|1 shows, the irreducible mass is a function of both the mass and the spin, and 
in general we have no completely unambiguous way specify the spin of the black holes in 
interaction. But, as was shown in 2y|, the spin evolves only very little until the black 
holes are very close together. Further, the apparent horizon coincides closely with the event 
horizon unless the black holes have strong interaction. Hence we assume that the individual 
spins are correctly given by the spin parameters ^^a specified in forming the superposed 
Kerr-Schild background, and that the mass is that determined by Eq. ()3()|1 using the area 
determined from the apparent horizon area. 

The physical idea in determining the binding energy is that the configuration is assembled 
from infinitely separated black holes, which are initially on the x-axis and which initially 
have parallel spins. (No energy is required to orient the coordinate system or to adiabatically 
rotate the spins while the holes are infinitely separated.) Thus these separated holes have 
their isolated total energy, i.e. Im^ for equal mass black holes. 

Then one of the black holes is adiabatically brought to a particular distance from the 
other, for instance a coordinate distance of lOrra as in some of our examples. This is the 
base configuration from which our computations will start. We then consider the change in 
the binding energy as we move the direction of the spin axis of one of the black holes. 

[|:] Dain [36j considers black hole slicings that have a second asymptotically flat infinity, and measures a 
mass (an intrinsic mass for the black hole) at this second infinity. This approach is impossible for the 
Kerr-Schild data we consider because Kerr-Schild slices intersect the black hole singularity. 
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9.2. ADM angular momentum 

Besides the mass, ADM formulae also exist for the momentum p^^^ and angular momen- 
tum JaiP^- These formulae are also evaluated in an asymptotically flat region surrounding 
the system of interest J33|, |37| : 



P 



ADM 



^j{Kki-K\5ki)dS\ (31) 

In the data we set, the total momentum is set to zero, so p^^^ = 0. In general, we set 
data for arbitrarily spinning holes with arbitrary orbital impact parameter, so in general 
the angular momentum JaiP^ is nonzero, and interesting. In the results presented below as 
code tests, we seek initially non-moving black holes, so the total angular momentum J^b^^ 
is simply the sum of the intrinsic angular momenta, S^maO. 



9.3. Computational Results 

We carried out several series of computational experiments to investigate the spin-spin 
interaction. In particular we considered instantaneously nonmoving black holes of equal mass 
rrii = 7712 = f^, with equal spin parameter ai = a2 = 0.5m. The background separation 
d for each series was varied from 6?7i to 18rra. For instance, we considered d = IOtti (holes 
at coordinate location xi = —5m, X2 = 5m). We varied the spin axis of one hole in two 
different planes, resulting in two "series" of data. The hole at x = +5m ( "hole number 2" ) 
was maintained with spin 02 aligned with the 2;-axis, while the direction of the other spin 
ai was varied in a plane in 7r/8 steps through 2tt from the -l-z-axis through the — 2;-axis and 
on back to the +2;-axis. The difference in the two series is that in one case (the "1/2; series" ) 
the spin remains in the y-z plane; in the other (the ^^xz series") it remains in the x-z plane. 
These two configurations are displayed in Figure 

The domains we used were typically ±15??7,, using 513'^ grid points, and typically excising 
a region of size rgx = 0.9??7,. The ADM mass was evaluated on a cube with sides at ±12m 
{i.e., well inside "convergence region" shown in Figured). (Variations in domain size, reso- 
lution, and excision region size were conducted to estimate the dependence of the resulting 
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FIG. 5: The two different BBH configurations investigated. In all cases, the black hole at x = +d/2 
is held fixed with a constant spin of 02 = 0.5 in the z direction. In the "xz series," the spin axis 
of hole at X = — d/2 is rotated in the x-z plane (i.e. about the y axis) by varying angle 0\ away 
from the other hole (holding ipi = it/2). In the "yz series," this axis is rotated in the y-z plane 
by varying 9i clockwise about the positive x axis. We note that, as an historical artifact of the 
background generator code of 2a], angles are defined such that y? = corresponds to the y axis. 



not the (more typical) x axis. 

binding energy on these physically irrelevant but computationally important parameters. 
For example we conducted a series of runs with outer boundaries at ±20?7i with 513^ points, 
evaluating the ADM mass at ±17m.) The apparent horizon areas were determined using 
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computational toolkit, via 



Thornburg's horizon finder [5J] in the Cactus 
a post-processing run on our output files. 

As shown in Figures |B] and [7| below, the angular dependence for the binding energy 
behavior in the yz case is close to that predicted by Wald (Eq. (|2ZI))- 

Figure El contains two tests of the yz series, computed identically except for a change in 
the excision radius. We see that the spin dependence of the binding energy is unchanged, 
but there is an offset in the average binding energy. This binding energy offset (0.001m out 
of — 0.07?7i) is a well known - but small - dependence on the inner boundary condition in the 
computation of initial data sets for binary black hole systems (see, e.g. 18|, J26|). It implies 
an accuracy in the binding energy (estimated from the value of the offset) of less than 2% 
of the total binding energy. On the other hand, the behavior of the spin-dependence, the 
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FIG. 6: (Normalized) Binding energy vs. spin angle for the yz series. Present in the graph are 
two curves corresponding to different excision radii Vex- For instance a least-squares fit to the 
Tex = 0.9m curve is Eb/{MAm + Mah2) = -0.06794 - 1.396 x 10-4cos6'. For the r^x = 0.75m 
curve, the amplitude of the cosine is 1.381 x 10~^. This cosine corresponds to the S ■ S' term in 
(|27|) . We note that changing the excision radius changes the overall constant offset of the binding 
energy, but does not have a large effect on the amplitude of the spin-spin interaction. 

cosine curve in the binding energy, implies a precision much smaller than the peak-peak 
amplitude of the cosine curve; we estimate 0.02%, one tenth of the peak to peak amplitude. 
However, Figure [3 for the xz series (where one of the spins tips away from the other) 
reminds us that there are additional physical effects in play. The Kerr solution has a 
quadrupole moment arising nonlinearly from its spin. In terms of Newtonian potential for 
Kerr: 



'^ ^ 2 2r^ 

'—r ~ iTT^^f^ cos B + ...; 
a 2# 



(33) 



the quadrupole term is the cos^ 6 term, where G is the "viewing" angle at which one hole 
"sees" the other. Given our configurations in which the spins by default are perpedicular to 
the line of separation, 6 = 7r/2 when ^ = 0. Thus in terms of our spin-orientation angle 6 
the effect may be written as 

m 3 



d 2d^ 



ma sin 9 + ...; 



(34) 
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FIG. 7: (Normalized) Binding energy vs. spin angle for the xz series ('Vi = '^/'^")j along with 
the Tex = 0.9m curve from Figure shown for comparison. Also included is a least-squares fit to 
the xz series: Ei,/{Mahi + Mah2) = -0.0679 - 1.390 x IQ-"^ cos 6 - 2.762 x 10""^ sin^ e. Note that 
the coefficient of the sin^ 9 term is roughly twice that of those cos 9 terms both in this graph an 
Figure El Rather than this near factor of 2 we find numerically, considerations based on the mass 
quadrupole moment would suggest a factor closer to 3/2, as indicated by (|51)) . 



Here d is a radial coordinate, defined so that angular dependence begins in the metric only 
at O((i~^)J30|. Hernandez [31| expands the asymptotic Kerr-Schild form and comes to the 
same result for the quadrupole moment of the Kerr black hole. The quadrupole cos^ G effect 
is not evident in the yz series because the hole at +5?Ti is always "looking" at the equator 
of the hole at — 5?7i, i.e. at G = 7r/2 so there is zero effect, even as the hole at —5m tilts. 
However, since the xz series tilts the hole at x = —5m away from hole at a; = +5m, the 
fixed hole "sees" different latitudes of the rotated hole in the xz series. 

Bonning et al. 26[ showed that Kerr-Schild data correctly predicts the Newtonian binding 
energy —mm'/d. The total binding in a relativistic calculation is this 0{d~^) term, plus 
Wald's 0{d~^) spin- spin interaction, plus the quadruole terms in the potential, plus any 
possible 0(o?~^) contribution to the solution. 

Following Wald's notation then, the complete spin dependence may be written as 
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Si ■ S2 — 3{Si ■ n){S2 ■ n)\ 3 /'"^2r^ -^-12 "^iro -12 



(Note that, for the configurations considered in this paper, S2 ■ n = 0.) Both the Wald 
spin-spin term and the quadrupole moment term in the expansion of the potential are 
proportional to d~^, though this is correct only near infinity; for close distances (such as 
d ~ lOm considered here) one expects deviations from nonlinear terms in the results. In fact 
the angular dependence of these terms is remarkably accurately reproduced. Table |l] shows 
the first /onr coefficients in fits to the binding energy; the third and fourth power coefficients 
are substantially below the cosine and cosine-squared coefficients. The Wald formula would 
produce an amplitude of 0.5^/10'^ = 2.5 x 10~^ for our case of equal spins of a = 0.5?Ti and 
separation oi d = 10m; the actual coefficient from the fit (after multiplying by the sum of 
the horizon masses) is 2.97 x 10""^. This apparent agreement is somewhat of an accident, 
however, since the expected dependence of d~^ is not present in our data, as we will show 
below. The term arising from the the quadrupole term (the cosine squared term) suggests 
a coefficient of 3.75 x 10~^ (1.5 times the expected amplitude of the spin-spin term). Our 
fit to the experiment (the xz series) produces 5.860 x 10"^. 

The Wald formula. Equation ()27|) . predicts no difference in the cosine term [Ai in Table 
H} between the xz and yz series. (The {S ■ n) term in Eq. (j^Tj) is zero for all experiments 
carried out because the black hole on the positive x axis has a fixed spin direction parallel 
to the z axis.) This is the behavior we find; compare the coefficients Bi for cos 6' in Tables 

mi and Em 

We tested the spin-squared dependence of the sin^ 6 term by two methods. In one we 
considered a = ai = 0.25?7i for the hole at x = —5m, which was then tested in an abbreviated 
yz series, while 02 held at 02 = 0.5 along the positive z axis for the hole at x = +5m. Figure 
IHl shows the result; the effect from the quadrupole term is quadratic in the reduced spin 
(its amplitude is reduced by a factor of four), while the Wald spin-spin interaction is linear 
in the reduced spin and its amplitude is reduced by a factor of two. To further test the 
quadrupole dependence we considered rotating the spin of the black hole at x = —5m in a 
plane turned by ipi = 0.175rad ~ 10°. The coefficients of the nonlinear curve fit are listed 
on the second line of Table H] and are plotted as interpolating lines in Figure [7| They have 
the analytically expected dependence. 
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Vl 
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^0 


Ai 


A2 


^3 


^4 


("yz series") 


10 


-0.0679461 


-1.3985e-4 


-1.07477e-6 


3.47174e-7 


-1.68858e-7 


0.175 


10 


-0.067933 


-1.3891e-4 


-1.45311e-5 


-1.39774e-6 


4.15175e-7 


1.57 ("xz series") 


10 


-0.0676704 


-1.3787e-4 


-2.74078e-4 


-1.55601e-6 


-2.17319e-6 



TABLE I: Table of coefficients for curve fits of the form Eb/{MAHi + Mah2) = ^o + ^i cos6' -H 
A2 cos^ 6 + A^ cos"^ 9 + A4 cos^ 9, for a separation of IOtti. (Here and below, we report several 
significant digits for the purposes of comparison, however due to variations resulting from excision 
region size and other effects, one would rightly regard only the first two digits as significant.) 
This shows that terms higher in order than cos^ 6 do not contribute significantly. Because of this, 
we do not include powers higher than two in the trigonometric basis functions. Also, given the 
considerations due to the mass quadrupole moment in Eq. ()34|) . all subsequent curve fits in this 
paper use the form Eh/{MAHi + Mah2) = Bq + Bi cos{9) + B2 sin 9. That is, the second order 
term will be taken as proportional to sin 9, not cos^ 9. This results in an offset of the total binding 
energy {Aq). 
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-0.000275675 


3.11827e-6 
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yz 
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-0.0576223 
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-8.578e-8 


yz 


16 


-0.0517695 


-6.45692e-5 


-3.17302e-7 


yz 


18 


-0.0459505 


-5.52391e-5 


-1.16e-6 



TABLE IL Table of coefficients for curve fits of the form Ei,/{Mahi + Mah2) = Bq + Bi cos9 + 
B2 SVC? 9, for the yz series {(pi = 0), for various BBH separations d with both spins ai = 02 = 0.5, 
using a domain size of ibl5m and 513^ fine grid points. Notice that, as expected, B2 is very small 
for all these fits. 
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FIG. 8: The effects of varying the spin magnitude of one of the holes. Symbols denote data points, 
lines denote curve fits. Notably, the fit for the xz series with ai = 0.5 is Ei,/{Mahi + Mah2) = 
-0.0679464 - 1.39041 x lO""^ cos 9 + 2.76249 x 10""^ sin^ 9, while for the xz series with m = 0.25 it 
is Eb/{MAHi + Mah2) = -0.0677621 - 6.94328 x 10"^ cos 9 + 7.00104 x 10^^ sin^ 9. Thus reducing 
spin ai from 0.5 to 0.25 results in reduction of the cos 9 term by a factor of two, while the sin^ 9 
term is reduced by nearly a factor of four. 
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TABLE III: Table of coefficients for curve fits of the form Ei,/{Mahi + Mah2) = Bq + Bi cos 9 + 
B2S\T?9, for the xz series ((/9i = 1.57), for various BBH separations d, using a domain size of 
iblSAf and 513^ fine grid points. 
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Figures IHl and ^1 show our tests of the separation-dependence of the binding energy. We 
expect the constant term Bq to fall off asymptotically as 1/d, since it corresponds to the 
M/r term in the Newtonian limit. Instead we find roughly linear behavior at the largest 
separations we are able to compute. The amplitudes Bi and B2 also fall off differently than 
expected. We expect both Bi, which corresponds to the cosine term in (j27|) . and B2, which 
corresponds to the mass quadrupole term, to scale as d~^. Instead we find that Bi scales 
as (i~^, and that B2 scales no faster than d~^. Since we expect the constant term Bq to 
scale as 1/d (although, as in Figure IHl we see that it does not), dividing the amplitudes Bi 
and B2 by the constant Bq does not significantly illuminate the results. These results for 
the separation-dependence are likely affected by the outer boundaries of our computational 
domain in unphysical ways. In the future we hope to repeat these studies with higher 
resolution and larger domains, using a multi-resolution (mesh refinement) version of our 
code. For the present, we conducted an additional test to measure the effects of the outer 
boundary, namely we looked for unphysical effects by computing the difference between the 
ADM mass and the horizon mass for a single black hole as we rotated its spin axis. The 
variation we found was on the order of 10~^, which would appear as a horizontal line in 
Figure El While this provides some assurance that our mass determination methods are 
functioning to some level of expectation, this (single-hole) effect is insufficient to explain 
the deviation from the expected results in the binary simulations. Future refinements and 
larger domains will hopefully clarify this issue. 

10. DISCUSSION 

To the extent checked, our computational tests of the spin-spin interaction in the binding 
energv of binary black hole configurations, verify the angular dependence given by Equation 
(I27|l p]. We verified the behavior given by Wald's —{S ■ S') term, but since we always 
kept one black hole's spin axis perpendicular to the separation axis, we made no attempt 
to observe the S ■ h term in Wald's formula. At the separations and domains available, our 
results did not show asymptotic d~^ fall-off with distance. 

Additionally, we find an effect due to the mass quadrupole moment of the black holes. 
This results in a higher-order (sine squared) variation with spin angle than given in Wald's 
formula. Thus we combine these two effects into a general equation (jHSI) for spin-spin 
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FIG. 9: Variation of binding energy vs. BBH separation d (holes at x = ±d/2), showing the 
constant term Bq in the curve fits shown in Table ^ The circles are computed using 513'^ grid 
points on a domain of ±15/71, evaluating the ADM mass on a cube of half-width 12m. We note that, 
at large separations, the curve is roughly linear, in contrast to an expected behavior of 1/d. We 
speculate that this (unphysical) effect is due to the outer boundary of the computational domain, 
in particular the behavior of the ADM mass as estimated at these rather small outer radii. However 
it may simply be the case that the asymptotic behavior only evident at larger separations. We 
expect that future work using our multiresolution code with larger domains and higher resolutions 
should produce better agreement with the expected behavior. 

interactions of binary black holes. However, in this case as well, the expected asymptotic 
fall off of d~^ was not evident in our solutions. Future analysis will use a new multiresolution 
version of our code to further pursue questions raised by the results here, including moving 
to significantly larger separations to evaluate the asymptotic behavior of the spin- dependent 
interactions, and effects due to rotation of both holes' spin axes. We also intend to investigate 
the spin-orbit coupling and its bearing on evolutions such as j^. We are now beginning 
exploration of the constrained evolution approach in spacetimes involving single moving, and 
multiple interacting black holes. We find substantial improvement from constraint solving 
in every simulation. 
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FIG. 10: Amplitude vs. separation for the constants Bi in Tables Hll and ITTll ('circles') . and B2 in 
Table UTII f squares). Rather than seeing the expected asymptotic d"^ fall-off for each amplitude, 
we find the cosine term Bi has roughly d~'^ behavior, whereas the sin^ term B2 falls off no faster 
than 1/d. 
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